/*
 * File:   mcFrustum.cpp
 * Author: Thales Luis
 *
 * Created on September 30, 2009, 5:33 PM
 */

#include "mcFrustum.h"

using namespace mc;

Frustum::Frustum() {
    // Inicializacao do quaternion de orientation
    this->orientation.set_value( 0.0f, 0.0f, 0.0f, 1.0f );

    // Inicializacao da perspectiva (Calcula os parametros)
    this->setPerspective( 30.0f, 1.0f, 0.0001f, 5000.0f );
}

Frustum::~Frustum() {}

// Posicao do observador: codifica a nova posicao em um quaternion com distancia
void Frustum::setPosition(float x, float y, float z) {
    // Atualiza posicao e distancia a origem
    this->x = x;
    this->y = y;
    this->z = z;

    this->computeParameters();
}

// Move o observador, dado um vetor
void Frustum::movePosition(float dx, float dy, float dz) {
      this->setPosition(x + dx, y + dy, z + dz);
}

// Avanca a posicao do observador
void Frustum::advancePosition(float frente, float acima, float direita) {
    this->setPosition(
        x + -view[0]*frente + right[0]*direita + up[0]*acima,
        y + -view[1]*frente + right[1]*direita + up[1]*acima,
        z + -view[2]*frente + right[2]*direita + up[2]*acima
    );
}


void Frustum::resetOrientation() {
    this->orientation.set_value( 0.0f, 0.0f, 0.0f, 1.0f );
    this->computeParameters();
}

// Configurador da direcao de visualizacao por vetor
void Frustum::setViewVector(float dirx, float diry, float dirz) {
    vec3<float> v(dirx, diry, dirz);
    vec3<float> dif;

    v = normalize(v);
    dif = v - this->view;

    if(length(dif) > EPSILON) {
        float  theta;
        quaternion<float> quat;
        vec3<float> n;

        // Find normal and angle cosine
        n.cross( this->view, v );
        theta = dot(v, this->view);

        // Check normal
        if( length(n) < 4 * EPSILON ) {
            n = this->right;
        } else {
            n = normalize(n);
        }

        // Check angle        
        if(theta > 1.0) theta = 1.0;
        else if(theta < -1.0) theta = -1.0;

        theta = acos(theta);

        /// Gira camera
        
        /// Axis to quaternion
        quat.set_value( n, theta );
        quat.normalize_quaternion();
        this->orientation *= quat;
        this->orientation.normalize_quaternion();

        this->computeParameters();
        this->view = v; // Force to avoid errors
    }
}

// Configurador da direcao de visualizacao com coordenadas geograficas (polares)
void Frustum::setViewVectorSpherical(float latit, float longit) {
    vec3<float> axis;
    quaternion<float> quat;

    latit  = deg_2_rad(latit);
    longit = deg_2_rad(longit);

    latit  = fmod((float) latit,  (float) M_PI);
    longit = fmod((float) longit, (float) M_PI);

    /// Reinicializa quaternion de orientacao
    this->orientation.set_value( 0.0f, 0.0f, 0.0f, 1.0f );

    /// Gira em Y
    axis[0] = 0.0f;
    axis[1] = 1.0f;
    axis[2] = 0.0f;    

    this->orientation *= quat;

    /// Gira em X
    axis[0] = 1.0f;
    axis[1] = 0.0f;
    axis[2] = 0.0f;

    /// Axis to quaternion
    quat.set_value(axis, -latit);
    this->orientation *= quat;
    this->computeParameters();
}


// Configurador do vetor cima por vetor
void Frustum::setUpVector(float up_x, float up_y, float up_z) {
    vec3<float> v(up_x, up_y, up_z );
    vec3<float> dif;

    v = normalize(v);

    dif = v - this->up;

    if( length(dif) > EPSILON ) {
        vec3<float> n;
        quaternion<float> quat;
        float theta;

        n.cross( this->up, v );
        theta = dot( v, this->up );

        if( length(n) < 4 * EPSILON )
            n = this->view;
        else
            n = normalize(n);

        if( theta > 1.0 )
            theta = 1.0;
        else if( theta < -1.0 )
            theta = -1.0;

        theta = acos(theta);

        quat.set_value( n, theta );
        quat.normalize_quaternion();
        this->orientation *= quat;
        this->orientation.normalize_quaternion();

        this->computeParameters();
        this->up = v;
    }
}

void Frustum::simulateTrackBall(
    quaternion<float> &quat,
    float p1x,
    float p1y,
    float p2x,
    float p2y
) {
    vec3<float> axis; /// Eixo de rotação
    float phi; /// Ângulo de rotação
    vec3<float> p1, p2, d;
    float t;

    if( p1x == p2x && p1y == p2y ) {
        quat.set_value( 0.0f, 0.0f, 0.0f, 1.0f );
        return;
    }

    p1.set_value(
        p1x,
        p1y,
        this->tb_project_to_sphere( 1.0, p1x, p1y )
    );

    p2.set_value(
        p2x,
        p2y,
        this->tb_project_to_sphere( 1.0, p1x, p1y )
    );

    axis.cross( p2, p1 );

    d = p1 - p2;
    t = length(d) / (2);

    if( t > 1.0 ) t = 1.0;
    if( t < -1.0 ) t = -1.0;
    phi = (float) asin(t);
    phi += phi;

    axis = normalize(axis);
    quat.set_value( axis, phi );
}



/**
 * Código de Gavin Bell
 * Projeta um par (x,y) em uma esfera de raio r OR uma folha hiperbólica
 * se estivermos longe do centro da esfera.
 */
inline float Frustum::tb_project_to_sphere( float r, float x, float y ) {
    float d, t, z;
    d = (float) sqrt(x*x + y*y);
   
    if( d < (float) (r * 0.70710678118654752440) )   /* Inside sphere */
        z = (float) sqrt(r*r - d*d);
    else {   /* On hyperbolic sheet */
        t = (float) (r / 1.41421356237309504880);
        z = t*t / d;
    }

    return z;
}





// Ajusta as direcoes para visualizar o alvo.
void Frustum::setTarget(float alvoX, float alvoY, float alvoZ) {
    this->setViewVector(x - alvoX, y - alvoY, z - alvoZ);
}


  // Funcao para girar a esfera que emcolve o visualizador
// Trackball: p1,p2 em [-1,1]x[-1,1] plano no qual a bola é projetada. rotacao livre.
void Frustum::rotateOrientationTrackball(float x1, float y1, float x2, float y2) {
    quaternion<float> quat; // Vetor velocidade angular da bola em forma de quartenion

    this->simulateTrackBall(quat, x1, y1, x2, y2);
  
    /// Verificar se multiplicação de quaternions é comutativa
    this->orientation = quat * this->orientation;

    /// Atualiza vetores de direcao
    this->computeParameters();
}

// Gira o planetário usando deslocamentos em angulos geograficos.
// Rotacao livre em graus do planetario.
void Frustum::rotateOrientationSpherical(float dlat, float dlong) {
    vec3<float> axis;
    quaternion<float> quat;

    dlat  = deg_2_rad(dlat);
    dlong = deg_2_rad(dlong);

    /// Gira em X
    axis.set_value( MC_X_AXIS );
    quat.set_value( axis, -dlat );
    this->orientation = quat * this->orientation;

    /// Gira em y
    axis.set_value( MC_Y_AXIS );
    quat.set_value( axis, -dlat );
    this->orientation = quat * this->orientation;

    /// Gira em z
    axis.set_value( MC_Z_AXIS );
    quat.set_value( axis, -dlat );
    this->orientation = quat * this->orientation;

    this->computeParameters();
}

// Rotation with Euler angles using the LOCAL COORDINATE SYSTEM.
void Frustum::rotateOrientationEuler(float pitch, float yaw, float roll) {
    quaternion<float> qx, qy, qz, quat;

    pitch = deg_2_rad( pitch );
    yaw   = deg_2_rad( yaw );
    roll  = deg_2_rad( roll );

    qx.set_value( this->right, pitch );
    qy.set_value( this->up,    yaw   );
    qz.set_value( this->view,  roll  );

    quat = qy * qz;
    quat = qz * quat;
    this->orientation *= quat;

    this->computeParameters();
}

// Rotation of the planetarium around a given axis.
void Frustum::rotateOrientationAxis(float theta, float axisX, float axisY, float axisZ) {
    quaternion<float> quat;
    vec3<float> n;

    theta = deg_2_rad( theta );

    // Gira camera
    n.set_value( axisX, axisY, axisZ );
    n = normalize(n);
    quat.set_value( n, theta );
    quat.normalize_quaternion();
    this->orientation *= quat;
    this->orientation.normalize_quaternion();

    this->computeParameters();
}

// Alinha o cima com Y, e lado com X
void Frustum::alignAxis() {
    vec3<float> axis;
    quaternion<float> quat;
    float lat, lon;

    // Calcula latitude
    if( this->view[1] >  1.0 )
        this->view[1] = 1.0;

    if( this->view[1] < -1.0 )
        this->view[1] = -1.0;

    lat = (float) ( (M_PI_2) - acos(this->view[1]) ); // Latitude atual

    if( lat > (M_PI_2 - 0.0076) )
        lat = (float) ( M_PI_2 - 0.0076 ); // Simply a small number
    else if( lat < -(M_PI_2 - 0.0076) )
        lat = (float) - (M_PI_2 - 0.0076 );

    // Temos que testar os quadrantes para a longitude em [0, 2PI)
    if( fabs(this->view[0]) < EPSILON)
        this->view[0] = ( ( this->view[0] < 0.0 ) ? - EPSILON : EPSILON );

    lon = (float) (atan((double) view[2] / (double) view[0])); // Longitude atual

    if(view[0] < 0.0)
        lon += (float) M_PI_2;
    
    lon = (float) (M_PI_2 - lon); // Longitude atual

    // Gira em Y
    axis.set_value( MC_Y_AXIS );
    this->orientation.set_value( axis, lon );

    // Gira em X
    axis.set_value( MC_X_AXIS );
    quat.set_value( axis, -lat );
    this->orientation = quat * this->orientation;

    this->computeParameters();
}


// Gira todo o sistema em torno de uma esfera com centro determinado.
// Trackball: p1,p2 em [-1,1]x[-1,1] plano no qual a esfera é projetada
void Frustum::rotateOrbitTrackball(float centerX, float centerY, float centerZ, float x1, float y1, float x2, float y2) {
    quaternion<float> quat, q;
    vec3<float> r( x - centerX, y - centerY, z - centerZ );

    this->simulateTrackBall(quat, x1, y1, x2, y2);
    quat *= this->orientation;

    q = this->orientation.conjugate();
    quat = q * quat; // Desalinha o sistema da camera

    r = quat.apply_to_vec3(r);

    this->setPosition( r.x + centerX, r.y + centerY, r.z + centerZ );
}

// Gira livremente a bola em graus.
void Frustum::rotateOrbitSpherical(float centerX, float centerY, float centerZ, float dlat, float dlong) { // Gira livremente a bola em graus.
    quaternion<float> quat, rot;
    vec3<float> axis;
    vec3<float> r( x - centerX, y - centerY, z - centerZ );

    dlat  = deg_2_rad( dlat  );
    dlong = deg_2_rad( dlong );

    // Gira em X
    axis.set_value( MC_X_AXIS );
    rot.set_value( axis, -dlat ); // Positivo sobe: mao direita

    // Gira em  Y
    axis.set_value( MC_Y_AXIS );
    quat.set_value( axis, dlong );
    rot = quat * rot;

    rot *= this->orientation;
    quat = this->orientation.conjugate(); // Alinha c/ sistema da camera+rotacao
    rot = quat * rot; // Desalinha o sistema da camera

    r = rot.apply_to_vec3(r); // Rotaciona o vetor r
    
    this->setPosition( r.x + centerX, r.y + centerY, r.z + centerZ);
}

// Gira livremente a bola em graus.
void Frustum::rotateOrbitEuler(float centerX, float centerY, float centerZ, float pitch, float yaw, float roll) {
    quaternion<float> rot, quat;
    vec3<float> r( x - centerX, y - centerY, z - centerZ );

    rot.set_value( pitch, yaw, roll );

    rot *= this->orientation;
    quat = this->orientation.conjugate();
    rot *= quat * rot;

    r = rot.apply_to_vec3(r);
    this->setPosition( r.x + centerX, r.y + centerY, r.z + centerZ );
}


// Rotacao do sistema de observacao em torno de um eixo que passa por O e tem direcao dir
void Frustum::rotateAxis(float angle, float ox, float oy, float oz, float dirX, float dirY, float dirZ) {
    vec3<float> axis;
    vec3<float> r( x - ox, y - oy, z - oz );
    quaternion<float> quat;

    angle = deg_2_rad(angle);

    axis.set_value( dirX, dirY, dirZ );
    axis = normalize(axis);
    quat.set_value( axis, angle );
  
    r = quat.apply_to_vec3(r); // Rotaciona o vetor r

    this->setPosition(r.x + ox, r.y + oy, r.z + oz);
}



// Configura a projecao perspectiva
// Calcula a matriz de projecao perspectiva diretamente e recalcula os dados do FRUSTUM.
void Frustum::setPerspective(float fovh, float aspect, float nearplane, float farplane) {
    // Calcula os outros planos
    float top    = nearplane *  tan( deg_2_rad( fovh / 2.0 ) );
    float bottom = -top;
    float right  = aspect * top;
    float left   = -right;

    // Monta a matriz de perspectiva OPENGL
    this->proj_matrix._00 = ((float) 2.0 * nearplane) / (right - left);
    this->proj_matrix._01 = 0;
    this->proj_matrix._02 = 0;
    this->proj_matrix._03 = 0;

    this->proj_matrix._10 = 0;
    this->proj_matrix._11 = ((float) 2.0 * nearplane) / (top - bottom);
    this->proj_matrix._12 = 0;
    this->proj_matrix._13 = 0;

    this->proj_matrix._20 = (right + left) / (right - left);
    this->proj_matrix._21 = (top + bottom) / (top - bottom);
    this->proj_matrix._22 = -(farplane + nearplane) / (farplane - nearplane);
    this->proj_matrix._23 = (float) -1.0;

    this->proj_matrix._30 = 0;
    this->proj_matrix._31 = 0;
    this->proj_matrix._32 = -((float) 2.0 * farplane * nearplane) / (farplane - nearplane);
    this->proj_matrix._33 = 0;

    this->computeParameters();
}


// Configura a projecao ortográfica
// Calcula a matriz de projecao ortográfica diretamente e recalcula os dados do FRUSTUM.
void Frustum::setOrthographic(float left, float right, float top, float bottom, float nearplane, float farplane) {
    // Monta a matriz ortográfica OPENGL
    this->proj_matrix._00 = (float) 2.0 / (right - left);
    this->proj_matrix._01 = 0;
    this->proj_matrix._02 = 0;
    this->proj_matrix._03 = 0;

    this->proj_matrix._10 = 0;
    this->proj_matrix._11 = (float) 2.0 / (top - bottom);
    this->proj_matrix._12 = 0;
    this->proj_matrix._13 = 0;

    this->proj_matrix._20 = 0;
    this->proj_matrix._21 = 0;
    this->proj_matrix._22 = (float) -2.0 / (farplane - nearplane);
    this->proj_matrix._23 = 0;

    this->proj_matrix._30 = -(right + left) / (right - left);
    this->proj_matrix._31 = -(top + bottom) / (top - bottom);;
    this->proj_matrix._32 = -(farplane + nearplane) / (farplane - nearplane);
    this->proj_matrix._33 = 1;

    this->computeParameters();
}


////////////////////////////////////////////////////////////////////////////////////
// Determina se um ponto estaŽ inserido no frustum.
////////////////////////////////////////////////////////////////////////////////////
point_pos Frustum::isPointVisible(float x, float y, float z ) {
	//   Na equacao de um plano (A*x + B*y + C*z + D = 0), [A B C] eŽ a normal do plano
	// e D eŽ a distancia do plano `a origem. Se um ponto (x y z) estaŽ sobre o plano,
	// a equacao eŽ verdadeira (=0). Se o resultado eŽ negativo, o ponto estaŽ atras do
	// plano. Se eŽ positivo, esta na frente do plano.
	//   Um ponto esta dentro do frustum, se ele estaŽ na frente (eq >= 0) de todos seus
	// planos.

	// Testa se o ponto estaŽ em frente ao plano anterior
	if(nearA * x + nearB * y + nearC * z < -nearD) return MC_FRUSTUM_OUTSIDE;

	// Testa se o ponto estaŽ em frente ao plano posterior
	if(farA * x + farB * y + farC * z  < -farD) return MC_FRUSTUM_OUTSIDE;

	// Testa se o ponto estaŽ em frente ao plano direito
	if(rightA * x + rightB * y + rightC * z < -rightD) return MC_FRUSTUM_OUTSIDE;

	// Testa se o ponto estaŽ em frente ao plano esquerdo
	if(leftA * x + leftB * y + leftC * z < -leftD) return MC_FRUSTUM_OUTSIDE;

	// Testa se o ponto estaŽ em frente ao plano inferior
	if(bottomA * x + bottomB * y + bottomC * z < -bottomD) return MC_FRUSTUM_OUTSIDE;

	// Testa se o ponto estaŽ em frente ao plano superior.
	if(topA * x + topB * y + topC * z < -topD) return MC_FRUSTUM_OUTSIDE;

	// O ponto esta atras de todos os planos
	return MC_FRUSTUM_INSIDE;
}


////////////////////////////////////////////////////////////////////////////////////
// Determina a visibilidade de uma esfera no frustum (TOTAL, PARCIAL ou NAOVISIVEL)
////////////////////////////////////////////////////////////////////////////////////
point_pos Frustum::isSphereVisible(float x, float y, float z, float raio) {
	// O centro da esfera pode estar fora do frustum mas a esfera nao. Assim.
	// deve-se levar o raio em consideracao. Assim, se uma equacao produziu 3, o centro
	// estaŽ a uma distancia -3 fora do frustum. Mas se o raio for 5, a esfera ainda
	// estaŽ dentro do frustum. Uma esfera esta fora do plano se sua distancia
	// eŽ maior que o raio.
	point_pos result = MC_FRUSTUM_INSIDE;

	// Testa a esfera em relacao ao plano proximo
	float distancia = nearA * x + nearB * y + nearC * z + nearD;
	if(distancia < -raio) return MC_FRUSTUM_OUTSIDE;
	else if(distancia <= raio) result = MC_FRUSTUM_PARTIAL;

	// Testa a esfera em relacao ao plano direito
	distancia = rightA * x + rightB * y + rightC * z + rightD;
	if(distancia < -raio) return MC_FRUSTUM_OUTSIDE;
	else if(distancia <= raio) result = MC_FRUSTUM_PARTIAL;

	// Testa a esfera em relacao ao plano esquerdo
	distancia = leftA * x + leftB * y + leftC * z + leftD;
	if(distancia < -raio) return MC_FRUSTUM_OUTSIDE;
	else if(distancia <= raio) result = MC_FRUSTUM_PARTIAL;

	// Testa a esfera em relacao ao plano inferior
	distancia = bottomA * x + bottomB * y + bottomC * z + bottomD;
	if(distancia < -raio) return MC_FRUSTUM_OUTSIDE;
	else if(distancia <= raio) result = MC_FRUSTUM_PARTIAL;

	// Testa a esfera em relacao ao plano superior
	distancia = topA * x + topB * y + topC * z + topD;
	if(distancia < -raio) return MC_FRUSTUM_OUTSIDE;
	else if(distancia <= raio) result = MC_FRUSTUM_PARTIAL;

        // Testa a esfera em relacao ao plano posterior (Menor probabilidade)
        distancia = farA * x + farB * y + farC * z + farD;
	if(distancia < -raio) return MC_FRUSTUM_OUTSIDE;
	else if(distancia <= raio) result = MC_FRUSTUM_PARTIAL;

	// A esfera está pelo menos parcialmente inserida no frustum
	return result;
}

////////////////////////////////////////////////////////////////////////////////////
// Determina a visibilidade de uma caixa alinhada com os eixos AABB.
////////////////////////////////////////////////////////////////////////////////////
point_pos Frustum::isBoxVisible(float xmin, float ymin, float zmin, float xmax, float ymax, float zmax) {
    // Usa o método de vértices negativos e positivos
    register float   xp, yp, zp;   // Coordenadas do vertice positivo
    register float   xn, yn, zn;   // Coordenadas do vertice negativo
    float   distance;

    ///////////////////////////////////////////////////////////////////////
    // Testa todos os vertices POSITIVOS (que devem estar dentro do plano)

    // Testa a caixa em relacao ao PLANO PROXIMO
    // Verifica a normal do plano e ajusta o vertice positivo (que deve estar dentro do plano)
    xp = (nearA >= 0) ? xmax : xmin;
    yp = (nearB >= 0) ? ymax : ymin;
    zp = (nearC >= 0) ? zmax : zmin;
    distance = nearA * xp + nearB * yp + nearC * zp + nearD;
    if(distance < 0) return MC_FRUSTUM_OUTSIDE; // Neste caso está totalmente fora do plano

    // Testa a caixa em relacao ao PLANO DIREITO
    // Verifica a normal do plano e ajusta o vertice positivo (que deve estar dentro do plano)
    xp = (rightA >= 0) ? xmax : xmin;
    yp = (rightB >= 0) ? ymax : ymin;
    zp = (rightC >= 0) ? zmax : zmin;
    distance = rightA * xp + rightB * yp + rightC * zp + rightD;
    if(distance < 0) return MC_FRUSTUM_OUTSIDE; // Neste caso está totalmente fora do plano

    // Testa a caixa em relacao ao PLANO ESQUERDO
    // Verifica a normal do plano e ajusta o vertice positivo (que deve estar dentro do plano)
    xp = (leftA >= 0) ? xmax : xmin;
    yp = (leftB >= 0) ? ymax : ymin;
    zp = (leftC >= 0) ? zmax : zmin;
    distance = leftA * xp + leftB * yp + leftC * zp + leftD;
    if(distance < 0) return MC_FRUSTUM_OUTSIDE; // Neste caso está totalmente fora do plano

    // Testa a caixa em relacao ao PLANO INFERIOR
    // Verifica a normal do plano e ajusta o vertice positivo (que deve estar dentro do plano)
    xp = (bottomA >= 0) ? xmax : xmin;
    yp = (bottomB >= 0) ? ymax : ymin;
    zp = (bottomC >= 0) ? zmax : zmin;
    distance = bottomA * xp + bottomB * yp + bottomC * zp + bottomD;
    if(distance < 0) return MC_FRUSTUM_OUTSIDE; // Neste caso está totalmente fora do plano

    // Testa a caixa em relacao ao PLANO SUPERIOR
    // Verifica a normal do plano e ajusta o vertice positivo (que deve estar dentro do plano)
    xp = (topA >= 0) ? xmax : xmin;
    yp = (topB >= 0) ? ymax : ymin;
    zp = (topC >= 0) ? zmax : zmin;
    distance = topA * xp + topB * yp + topC * zp + topD;
    if(distance < 0) return MC_FRUSTUM_OUTSIDE; // Neste caso está totalmente fora do plano


    // Testa a caixa em relacao ao PLANO DISTANTE
    // Verifica a normal do plano e ajusta o vertice positivo (que deve estar dentro do plano)
    xp = (farA >= 0) ? xmax : xmin;
    yp = (farB >= 0) ? ymax : ymin;
    zp = (farC >= 0) ? zmax : zmin;
    distance = farA * xp + farB * yp + farC * zp + farD;
    if(distance < 0) return MC_FRUSTUM_OUTSIDE; // Neste caso está totalmente fora do plano


    // Neste ponto sabemos que a caixa esta parcialmente ou completamente dentro do FRUSTUM

    ///////////////////////////////////////////////////////////////////////
    // Testa todos os vertices NEGATIVOS (que devem estar fora do plano)

    // Testa a caixa em relacao ao PLANO PROXIMO
    // Verifica a normal do plano e ajusta o vertice negativo (que deve estar fora do plano)
    xn = (nearA >= 0) ? xmin : xmax;
    yn = (nearB >= 0) ? ymin : ymax;
    zn = (nearC >= 0) ? zmin : zmax;
    distance = nearA * xn + nearB * yn + nearC * zn + nearD;
    if(distance < 0) return MC_FRUSTUM_PARTIAL; // Neste caso sabemos que é ao menos PARCIAL

    // Testa a caixa em relacao ao PLANO DIREITO
    // Verifica a normal do plano e ajusta o vertice negativo (que deve estar fora do plano)
    xn = (rightA >= 0) ? xmin : xmax;
    yn = (rightB >= 0) ? ymin : ymax;
    zn = (rightC >= 0) ? zmin : zmax;
    distance = rightA * xn + rightB * yn + rightC * zn + rightD;
    if(distance < 0) return MC_FRUSTUM_PARTIAL; // Neste caso sabemos que é ao menos PARCIAL

    // Testa a caixa em relacao ao PLANO ESQUERDO
    // Verifica a normal do plano e ajusta o vertice negativo (que deve estar fora do plano)
    xn = (leftA >= 0) ? xmin : xmax;
    yn = (leftB >= 0) ? ymin : ymax;
    zn = (leftC >= 0) ? zmin : zmax;
    distance = leftA * xn + leftB * yn + leftC * zn + leftD;
    if(distance < 0) return MC_FRUSTUM_PARTIAL; // Neste caso sabemos que é ao menos PARCIAL

    // Testa a caixa em relacao ao PLANO INFERIOR
    // Verifica a normal do plano e ajusta o vertice negativo (que deve estar fora do plano)
    xn = (bottomA >= 0) ? xmin : xmax;
    yn = (bottomB >= 0) ? ymin : ymax;
    zn = (bottomC >= 0) ? zmin : zmax;
    distance = bottomA * xn + bottomB * yn + bottomC * zn + bottomD;
        if(distance < 0) return MC_FRUSTUM_PARTIAL; // Neste caso sabemos que é ao menos PARCIAL

    // Testa a caixa em relacao ao PLANO SUPERIOR
    // Verifica a normal do plano e ajusta o vertice negativo (que deve estar fora do plano)
    xn = (topA >= 0) ? xmin : xmax;
    yn = (topB >= 0) ? ymin : ymax;
    zn = (topC >= 0) ? zmin : zmax;
    distance = topA * xn + topB * yn + topC * zn + topD;
    if(distance < 0) return MC_FRUSTUM_PARTIAL; // Neste caso sabemos que é ao menos PARCIAL

    // Testa a caixa em relacao ao PLANO DISTANTE
    // Verifica a normal do plano e ajusta o vertice negativo (que deve estar fora do plano)
    xn = (farA >= 0) ? xmin : xmax;
    yn = (farB >= 0) ? ymin : ymax;
    zn = (farC >= 0) ? zmin : zmax;
    distance = farA * xn + farB * yn + farC * zn + farD;
    if(distance < 0) return MC_FRUSTUM_PARTIAL; // Neste caso sabemos que é ao menos PARCIAL

    // A caixa está no frustum totalmente
    return MC_FRUSTUM_INSIDE;
}

// Usa compilação condição para evitar link com bibliotecas do OPENGL caso o usuário não queira usar
#ifdef OPENGL
    // Exporta as transformacoes para o OpenGL
    void Frustum::exportOpenGL() {
        glMatrixMode(GL_PROJECTION);
        #ifndef USE_GL_DOUBLE
            glLoadMatrixf(proj_matrix);
            glMatrixMode(GL_MODELVIEW);
            glLoadMatrixf(view_matrix);
        #else
            glLoadMatrixd(proj_matrix);
            glMatrixMode(GL_MODELVIEW);
            glLoadMatrixd(view_matrix);
        #endif
    }
#endif



// Calcula os vetores de visualizacao a partir do quaternion
// Calcula os planos do frustum
void Frustum::computeParameters() {
    matrix4<float> m, rot, trans;
    matrix4<float> transM; // Guarda a transformacao completa TRANSPOSTA(projecao * modelovisao)

    register float A, B, C, D, norma;


    // Extrai a orientacao da camera do quaternion
    this->orientation.normalize_quaternion();
    this->orientation.get_value(m); // Padrao OPENGL
    
    this->right.set_value( m._00, m._10, m._20 );
    this->up.set_value(    m._01, m._11, m._12 );
    this->view.set_value(  m._02, m._12, m._22 );


    // Trying to minimize error propagation!
    this->view  = normalize(this->view);
    this->up    = normalize(this->up);
    this->right = normalize(this->right);

    m._00 = this->right.x;
    m._01 = this->up.x;
    m._02 = this->view.x;

    m._10 = this->right.y;
    m._11 = this->up.y;
    m._12 = this->view.y;

    m._20 = this->right.z;
    m._21 = this->up.z;
    m._22 = this->view.z;

    rot.transpose(m);
  

    // Monta a matriz de translacao (TRANSPOSTA DO PADRAO OPENGL)
    trans._03 = -this->x;
    trans._13 = -this->y;
    trans._23 = -this->z;

    // Calcula a matriz de visualizacao
    m = rot * trans;
    m.transpose( this->view_matrix ); // Padrao OpenGL

    // Multiplica as transformacoes de projecao e modelo_visao que estao transpostas
    // (padrao OPENGL). De acordo com a propriedade de matrizes (A . B)t = Bt . At,
    // a transposta da transformacao M = (projecao . modelo_visao) eŽ obtido pela
    // multiplicacao de TRANSPOSTA(M) = (TRANSPOSTA(modelo_visao) . TRANSPOSTA(projecao)).

    transM = this->view_matrix * this->proj_matrix;

    //////////////////////////////////////////////////////////////////////////
    //	  Por default, uma camera o OPENGL eŽ direcionada para o z negativo com
    // o lado superior da cabeca voltada para o eixo y positivo. Assim, inicialmente
    // o vetores normais (coordenadas homogeneas) aos planos do volume sao:
    //		- plano superior  = [ 0 -1  0  1 ]
    //		- plano inferior  = [ 0  1  0  1 ]
    //		- plano posterior = [ 0  0  1  1 ]
    //		- plano anterior  = [ 0  0 -1  1 ]
    //		- plano direito   = [-1  0  0  1 ]
    //		- plano esquerdo  = [ 1  0  0  1 ]
    //
    //	  Esses vetores estao imcertidos para que a tranformacao de projecao calcule
    // corretamente os respectivos vetores normais dos planos que apontam para dentro
    // do frustum;
    //    Os planos do frustum sao obtidos pela tranformacao dos vetores normais
    // aos planos. Com isso obtem-se a equacao dos planos Ax + By + Xz + D = 0, onde
    // [A B C] eŽ o vetor normal ao plano e D eŽ a sua distancia `a origem.

    ///////////////////////////////////////////////////////////////
    // Obtem a equacao do plano superior do volume = M * [0 -1 0 1]
    ///////////////////////////////////////////////////////////////

    A = transM._03 - transM._01;
    B = transM._13 - transM._11;
    C = transM._23 - transM._21;
    D = transM._33 - transM._31;

    // Normaliza o plano superior dividindo a sua equacao pela norma do vetor [A B C]
    norma = (float) sqrt(A * A + B * B + C * C);
    this->topA = A / norma;
    this->topB = B / norma;
    this->topC = C / norma;
    this->topD = D / norma;

    ///////////////////////////////////////////////////////////////
    // Obtem a equacao do plano inferior do volume = M * [0 1 0 1]
    ///////////////////////////////////////////////////////////////
    A = transM._03 + transM._01;
    B = transM._13 + transM._11;
    C = transM._23 + transM._21;
    D = transM._33 + transM._31;

    // Normaliza o plano inferior dividindo a sua equacao pela norma do vetor [A B C]
    norma = (float) sqrt(A * A + B * B + C * C);
    this->bottomA = A / norma;
    this->bottomB = B / norma;
    this->bottomC = C / norma;
    this->bottomD = D / norma;

    ///////////////////////////////////////////////////////////////
    // Obtem a equacao do plano posterior do volume = M * [0 0 1 1]
    ///////////////////////////////////////////////////////////////
    A = transM._03 + transM._02;
    B = transM._13 + transM._12;
    C = transM._23 + transM._22;
    D = transM._33 + transM._32;
    

    // Normaliza o plano posterior dividindo a sua equacao pela norma do vetor [A B C]
    norma = (float) sqrt(A * A + B * B + C * C);
    this->farA = A / norma;
    this->farB = B / norma;
    this->farC = C / norma;
    this->farD = D / norma;

    ///////////////////////////////////////////////////////////////
    // Obtem a equacao do plano anterior do volume = M * [0 0 -1 1]
    ///////////////////////////////////////////////////////////////
    A = transM._03 + transM._02;
    B = transM._13 + transM._12;
    C = transM._23 + transM._22;
    D = transM._33 + transM._32;

	// Normaliza o plano anterior dividindo a sua equacao pela norma do vetor [A B C]
    norma = (float) sqrt(A * A + B * B + C * C);
    this->nearA = A / norma;
    this->nearB = B / norma;
    this->nearC = C / norma;
    this->nearD = D / norma;

    ///////////////////////////////////////////////////////////////
    // Obtem a equacao do plano direito do volume = M * [-1 0 0 1]
    ///////////////////////////////////////////////////////////////
    A = transM._03 - transM._00;
    B = transM._13 - transM._10;
    C = transM._23 - transM._20;
    D = transM._33 - transM._30;

    // Normaliza o plano direito dividindo a sua equacao pela norma do vetor [A B C]
    norma = (float) sqrt(A * A + B * B + C * C);
    this->rightA = A / norma;
    this->rightB = B / norma;
    this->rightC = C / norma;
    this->rightD = D / norma;

    ///////////////////////////////////////////////////////////////
    // Obtem a equacao do plano esquerdo do volume = M * [1 0 0 1]
    ///////////////////////////////////////////////////////////////
    A = transM._03 + transM._00;
    B = transM._13 + transM._10;
    C = transM._23 + transM._20;
    D = transM._33 + transM._30;


    // Normaliza o plano esquerdo dividindo a sua equacao pela norma do vetor [A B C]
    norma = (float) sqrt(A * A + B * B + C * C);
    this->leftA = A / norma;
    this->leftB = B / norma;
    this->leftC = C / norma;
    this->leftD = D / norma;
}